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Abstract 

Many time-dependent deformation processes at elevated temperatures produce significant concurrent mi¬ 
crostructure changes that can alter the mechanical properties in a profound manner. Such microstructure 
evolution is usually absent in mesoscale deformation models and simulations. Here we present an inte¬ 
grated full-held modeling scheme that couples the mechanical response with the underlying microstruc¬ 
ture evolution. As a hrst demonstration, we integrate a fast Fourier transform-based elasto-viscoplastic 
(FFT-EVP) model with a phase-held (PF) recrystallization model, and carry out three-dimensional sim¬ 
ulations of dynamic recrystallization (DRX) in polycrystalline copper. A physics-based coupling between 
FFT-EVP and PE is achieved by (1) adopting a dislocation-based constitutive model in EET-EVP, which 
allows the predicted dislocation density distribution to be converted to a stored energy distribution and 
passed to PE, and (2) implementing a stochastic nucleation model for DRX. Calibrated with the ex¬ 
perimental DRX stress-strain curves, the integrated model is able to deliver full-held mechanical and 
microstructural information, from which quantitative description and analysis of DRX can be achieved. 
It is suggested that the initiation of DRX occurs signihcantly earlier than previous predictions, due to 
heterogeneous deformation. DRX grains are revealed to form at both grain boundaries and junctions 
(e.g., quadruple junctions) and tend to grow in a wedge-like fashion to maintain a triple line (not neces¬ 
sarily in equilibrium) with old grains. The resulting stress redistribution due to strain compatibility is 
found to have a profound inhuence on the subsequent dislocation evolution and softening. 

Keywords: Dynamic recrystallization. Thermomechanical processes. Crystal plasticity. Phase-held, 
Microstructures 


1. Introduction 

Many thermomechanical processes and high temperature applications of materials involve a coupled 
evolution of micromechanical helds, local defect populations, and microstructural constituents including 
precipitates and grains. Eor instance, during hot deformation crystalline materials with low stacking 
fault energy (SEE) often undergo dynamic recrystallization (DRX) wherein new grains will continue to 
nucleate and grow (Sakai and Jonas, 1984; Sakai et ah, 2014), altering the population, mutual elastic in¬ 
teraction, and subsequent motion of dislocations in a distinctive manner. These microstructural changes 
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are difficult to capture through purely mechanical laws (Roters et ah, 2010). The continued adoption of 
new computational strategies for accelerated development of materials, such as Integrated Computational 
Materials Engineering (Allison et ah, 2006; Allison, 2011), for materials processed through thermome¬ 
chanical routes or materials exposed to thermal and mechanical extremes in-service requires modeling 
and simulation tools that integrate both the mechanical and microstructural aspects in a fully coupled 
manner. 

Within the broader solid mechanics and materials science communities, mature computational approaches 
exist to study, separately, the deformation of heterogenous materials and microstructure evolution, re¬ 
spectively. On one hand, the field of crystal plasticity has benefited substantially from multi-scale 
experiments and simulations linking mechanical properties of materials with the evolution of contained 
structural defects such as dislocations under an applied stress or displacement field. At the mesoscale, 
this understanding has been implemented into physics-based constitutive theories (Arsenlis and Parks, 
2002; Arsenlis et ah, 2004; Cheong and Busso, 2004; Ma et ah, 2006; Gao and Huang, 2003; Beyerlein 
and Tome, 2008) and implemented into homogenized deformation models such as self-consistent schemes 
(Lebensohn and Tome, 1993; Niezgoda et ah, 2014) or full-field simulations such as finite element based 
crystal plasticity (FE-CP) (Kalidindi et ah, 1992; Beaudoin et ah, 1995; Roters et ah, 2010) or fast Fourier 
transform (EFT) based crystal plasticity (FFT-CP) models (Lebensohn, 2001; Lebensohn et ah, 2012; 
Eisenlohr et ah, 2013). On the other hand, the microstructural evolution in crystals, such as grain growth 
(Chen and Yang, 1994; Kazaryan et ah, 2002; Moelans et ah, 2008b), static recrystallization (Moelans 
et ah, 2013), rafting in superalloy (Zhou et ah, 2010; Gaubert et ah, 2010) and many other phenomena 
(Chen, 2002; Wang and Li, 2010) have been well studied using phase-field (PE) simulations. The non¬ 
boundary tracking field description of microstructures and the incorporation of thermodynamics-based 
free energy formulation have made PE a very powerful and robust tool in simulating and predicting 
the microstructural evolution often in a quantitative manner (Chen, 2002; Boettinger et ah, 2002; Shen 
et ah, 2004; Moelans et ah, 2008b; Steinbach, 2009; Wang and Li, 2010). 

In recent years, efforts have been made towards developing CP models that can incorporate microstruc¬ 
ture features to study plastic deformation of materials such as austenitic steels, TRIP steels, brass, 
TWIP steels and shape memory alloys that deform not only by dislocation slips but also by displacive 
phase transformation mechanisms. For instance, frameworks (Thamburaja and Anand, 2001; Turteltaub 
and Suiker, 2005; Lan et ah, 2005; Manchiraju and Anderson, 2010) have been developed to incorporate 
martensitic transformations as the flow rules. Mechanical twinning, which is of great importance to the 
plasticity of many BCC metals as well as FCC metals with low SEE, has also been incorporated into 
FE-CP models (Kalidindi, 1998; Staroselsky and Anand, 1998; Salem et ah, 2005; Steinmetz et ah, 2013; 
Zhang et ah, 2008). Atomistically-informed dislocation-based models have also been developed recently 
and applied to single crystal (Cereceda et ah, 2015). In addition, microstructure modeling techniques 
such as cellular automata and phase-field have also been applied to the study of mechanics-induced 
microstructural evolution such as static recrystallization (Hesselbarth and Gobel, 1991; Raabe, 2002; 
Takaki et ah, 2007; Moelans et ah, 2013; Chen et ah, 2015). These models, while providing certain 
connection between microstructure and crystal plasticity, still lack a dynamic coupling between the two. 
On the other hand, phase-field models incorporating either continuum plasticity (Gaubert et ah, 2010) 
or dislocation density fields (Zhou et ah, 2010) have also been developed to study rafting in Ni-based 
superalloys. These models are mainly rooted in the PE framework and significant advances are required 
in order to generalize the approach and incorporate a wider range of existing constitutive theories. 

Regarding the simulation of DRX, there have been models aiming to couple deformation with microstruc¬ 
ture evolution. In the models of Ding and Guo (2001) and Takaki et al. (2008), the growth of DRX grains 
was modeled by cellular automaton and PE, respectively, and the phenomenological Kocks-Mecking 
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model (Mecking and Kocks, 1981) was used to model the evolution of average dislocation density, which 
in return influenced the mechanics through its square root relationship with the flow stress. Recently 
Takaki et ah (2014) extended their previous work by replacing the flow stress model with an elastic- 
plastic FE, intending to establish a full coupling model. However, the mechanical behavior was still 
considered in a macroscopic level in the sense that an average dislocation density was assumed for each 
grain which evolved according to the Kocks-Mecking model (Mecking and Kocks, 1981). DRX, on the 
other hand, occurs at a sub-grain level where grain boundary (GB) bulging (Ponge and Gottstein, 1998; 
Wusatowska-Sarnek et ah, 2002; Miura et ah, 2007) or other possible mechanisms (Rollett et ah, 2004) 
might operate to initiate the nucleation of new grains. This length scale separation implies that a “ho¬ 
mogenized” coupling scheme (Takaki et ah, 2014) may limit its application in exploring the underlying 
physics. Another recent DRX model by Popova et al. (2014), which couples FE-CP with a probabilistic 
cellular automata, successfully captured both texture and mechanical feature during the DRX in mag¬ 
nesium. However, only geometrically necessary dislocations (GND) were considered with the effect of 
statistically stored dislocations (SSD) being ignored, and the two-dimensional (2D) simulation ignores 
important microstructural features (e.g., the triple/quadruple grain junctions (Miura et ah, 2005)) when 
analyzing the dynamics of DRX grains. 

In this paper, we present an integrated modeling scheme of fully coupling the mechanical response 
with the underlying microstructure evolution by employing FFT-CP and PE. This novel scheme is then 
implemented in simulating dynamic recrystallization in polycrystalline copper in three-dimensions (3D). 
Quantitative agreement between simulations and experiments is achieved, and the revealed underlying 
DRX structures and evolution, to the best of our knowledge, are investigated in full-held 3D dynamics 
for the hrst time. 


2. Model 

To better describe the integrated model that fully couples the mechanics and microstructure, we need to 
consider both the abstract model structure (AMS), which describes how the PE and mechanical simulation 
environments interact, and a concrete model structure (GMS), which describes the specihc constitutive 
and kinetic models and what microstructural held variables are passed, as shown in Fig. 1. At the 
AMS level, the framework consists of three parts: (1) deformation kinematics, (2) description of model 
interfaces, including a constitutive theory that describes the evolution of a set of microstructural state 
variables (i.e., microscopic variables such as dislocation density and arrangement, cell size, grain size, 
precipitate size, and spacing, and so forth (Frost and Ashby, 1982)) and their interplay with PE order 
parameters, and (3) PE governing equations describing the kinetics of microstructure evolution as a 
function of order parameters. Additional nucleation models are needed when considering issues involving 
the formation of new grains (e.g., DRX and twinning) or new phases (e.g., precipitation) during the 
deformation. Implementation of the AMS requires the identihcation of microstructure descriptors as 
well as GP frameworks, constitutive laws, PE equations, and/or nucleation models, which all together 
describe the GMS. As our hrst demonstration, we have chosen modeling DRX in polycrystalline copper, 
of which the corresponding CMS with the basic simulation howchart is shown in Fig. 1(b). Based on 
the experiment of Ghauri et al. (1990), it was estimated that during a compression test up to 130% 
strain under applied strain rate of 1.6 x 10“^ /s (conditions to be used in the following simulation), grain 
growth in a polycrystal copper with an average grain size > 20 fim can be safely ignored (relative change 
< 2%). As a result, for efficiency and simplicity, in Fig. 1(b) PE is entered only when DRX is initiated. 
In the following, the details of each model element will be presented. 
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(a) (b) 


Figure 1: An integrated modeling scheme for thermal-mechanical processes: (a) a generic abstract model structure 
and (b) a concrete model structure showing the flowchart for simulating dynamic recrystallization. Note that 
once DRX is initiated, the nucleation occurs so frequently that the PF microstructure relaxation occurs virtually 
every FFT-EVP step (see more detailed discussion about implementing (b) in Sec. 2.5). 


2.1. Deformation kinematics 

We employ the fast Fourier transform-based elasto-viscoplastic (FFT-EVP) formulation in the framework 
of infinitesimal-strain, developed by Lebensohn et al. (2012). Using an Euler implicit time discretization 
and anisotropic elasticity, the stress in material point x at t + At is given by 

^t+At (x) = c(x) : ( 1 ) 

= C (x) ; [£*+^* (x) - £P-* - (x, At] 

where C is the elastic stiffness tensor, cr is the stress field, and the total strain field s, in the small-strain 
assumption, is simply the sum of elastic and plastic part, i.e. s = s® -f The plastic strain at t + At is 
obtained by the first-order approximation + ^p,t+At usage of the plastic 

strain rate at t + At (a function of stress at t + At) in calculating the plastic shear increment from t 
to t + At suggests an implicit Euler method, which requires the stress field to be solved iteratively, here 
by a spectral approach which utilizes FFT algorithms to simplify convolution integrals involving Green’s 
function. The plastic strain rate is formulated as 

AT 

£P(x) = ^m“(x) 7 “(x) (2) 

a=l 

where 7 ^ (x) and (x) are, respectively, the shear rate and the Schmid tensor associated with slip 
system a at point x, and Af is the number of active slip systems under consideration. The exact form of 
7 ^ is determined by the adopted constitutive laws (Sec. 2 . 2 ), which are usually nonlinear equations and 
will be solved with an augmented Lagrangian scheme adapted from Michel et al. (2000). 

The small-strain assumption in Eq. (1) implies that (i) the local neighborhood of a material point does 
not change drastically due to pure deformation, and (ii) the deformation increment is small with each 
numerical time step. These can be considered as legitimate in modeling DRX, of which the onset is still 
at a relatively small strain level (usually ^0.1) (Sakai and Jonas, 1984; Rollett et ah, 2004). In addition, 
once DRX is initiated, the mechanical behavior starts to be dominated by the migration of DRX grain 


4 


















boundaries modeled by PF, and the error due to small-strain assumption is further relieved. Additionally, 
the visco-plastic FFT (no elasticity) formulation with small-strain updates has been demonstrated by 
Prakash and Lebensohn (2009) to give comparable results to crystal-plasticity finite elements in the 
simulation of rolling in f.c.c. metals up to 60% thickness reduction and drawing in b.c.c. metals up to 40% 
equivalent strain. In that work it was found that even at large deformations the FFT scheme captured 
the location and distribution of stress fluctuations accurately, but that the specific numerical value of 
the stress localization varies with respect to the CP-FEM. These previous results provide confidence that 
stress/strain localizations predicted by the FFT-EVP will at least be statistically representative of those 
occurring in the material and sufficient to accurately predict evolution of the dislocation densities required 
for modeling nucleation and growth of new grains during DRX, even if the specific spatial values differ. 
Recently Eisenlohr et al. (2013) demonstrated a EET formulation which incorporates finite deformation 
kinematics, which utilizes the 1st Piola-Kirchhoff stress and the deformation gradient as work conjugate 
stress/deformation measures. Coupling of the finite strain EET with phase field equations poses several 
additional computational difficulties, namely that the Eourier grid for the finite strain EET is naturally 
formulated in the material or reference configuration while the spectral grid for the phase field is defined 
in the spatial or deformed state. Adoption of finite strain kinematics would require either a) explicitly 
mapping the material configuration to the spatial and interpolating onto a regular Eourier grid for the 
PE and performing the inverse mapping/interpolation every simulation time/strain increment or b) 
reformulating the phase field in the undeformed configuration. Another possible correction would be to 
adopt the approach of Cuitino and Ortiz (1992) which uses a logarithmic mapping to extend the small- 
strain updates to the finite deformation regime. However that approach requires the approximation that 
the plastic spin is identically zero and the principal directions of the plastic strain are aligned with the 
principal stresses. Eor these reasons we adopt the general EET-EVP framework of Lebensohn et al. 
(2012), keeping in mind the inherent limitations. Integration of the finite strain EET and phase field 
is ongoing and will be reported in future work. Specific discussion on the limitations and implications 
of the small strain kinematics with respect to the specific DRX case study investigated can be found in 
Sec. 3. 


2.2. Crystal plasticity constitutive laws 


The model interface should provide a methodology link between the deformation mechanics and mi¬ 
crostructure evolution such that the integrated modeling is physically self-consistent. This requires the 
shear rate in Eq. (2) to be formulated based on a set of state variables that can also be used as mi- 
crostructure descriptors. Eollowing the spirit of Orowan equation 7 = pu^v where pu is the mobile 
dislocation density, b the magnitude of Burgers vector, and v the average dislocation velocity, Ma and 
Roters (2004) and Ma et al. (2006) derived a dislocation-based flow rule for single-phase ECC metals: 
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where is the plastic shear rate for slip system a (among the total Af active slip systems) and the 
pre-exponential factor 70 = with the attempt frequency z/d being in the order of Debye 

frequency, /cb the Boltzmann constant, T the temperature, p the shear modulus, and q fitting parameters. 
Qsiip is the effective activation energy barrier for gliding through distributed obstacles consisting of 
parallel and forest dislocations with a density of pp and pp, respectively, which result in a passing stress 
^pass ~ and a cutting stress respectively. According to Ma and Roters (2004) 
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and Ma et al. (2006), pp and pp are given as 


Pv = £ 4 sD|cOs(ii“,F) I + Ip^nDsCOS | + |/)§NDet cos | 

/3=1 

+ IPGNDenCOS |, 


(4) 


/’P = £pSSDlsin(n“,t^) | + |pSNDsSin(n“,d^) | + |p§NDetSin(n“,t^) I 

/3=i 

+ I^GNDen sm (n“,n^) |, 


(5) 


where h^, and are, respectively, slip plane normal, slip direction and sense vector of slip system 
Pgnds’ PcNBet PcNDen related to the total GND density Pq^d system alpha and will 

be defined in a moment. The SSD density Pgsp evolves as (Ma et ah, 2006) 


PSSD = C4v/^ 7“ + C6ddipoiePM7“ “ C5 Pssd 7“ “ c? exp 


Qbulk 


knT 


kBT 


(PsSd) (7vm) 


( 6 ) 


which accounts for (1) the lock forming between mobile and forest dislocations, (2) the dipole formation 
between mobile dislocations, (3) the athermal annihilation of two parallel dislocations with anti-parallel 
Burgers vectors within a critical distance and (4) thermal annihilation by edge dislocation climb. Here 
G^dipoie = i 67 r^-t)r^ Critical distance for dipole formation (with u being the Poisson’s ratio), 

Qbuik the activation energy for self-diffusion, the von Mises equivalent stress and 7vm the von Mises 
equivalent shear rate. 


The GND density Pgnd? based on the crystallographic relationship, is decomposed into Pgnds’ PcNDet 
and PcNDen ^t ah, 2006). The rate of these effective GND densities satisfy the conservation law 
(Pgnd)^ = (Pgnds)^ + (PcNDet)^ + (PcNDen)^ cau be determined from the deformation gradient 
following the pioneering works done by Nye (1953), Dai and Parks (1997) and Dai (1997), which in small 
strain approximation reduce to (Dai, 1997) 


Pgnds = 

PGNDet = ■ d“, 

PCNDen = b* 


The gradient terms in Eq. (7) involve the physical size of simulation grid Iq, which is in principle 
determined by the physical length scale of the representative volume element (RVE). It is worth noting 
that when numerically evaluating these gradient terms, forward or backward difference will be used for 
grid points at grain boundaries and central difference for points in the grain interior. Galculation of GND 
through Eq. (7), however, ignore the effect that the grain boundary or interface, where the deformation 
gradient usually build, can also act as a source/sink to emit/absorb GND (Sun et ah, 2000). Recently, the 
saturation of GND has been studied both experimentally (Kysar et ah, 2010) and theoretically (Oztop 
et ah, 2013). In the current model, for simplicity, we account for this effect via a phenomenological 
treatment. In calculating the gradient in Eq. (7), instead of /q, we use an “effective” length /eff, updated 
as 

‘T‘ = 4 (i + x) ’ 

where ^vm is the current average von Mises equivalent strain and ^ (^ 1) is a fitting parameter that will 
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give rise to a saturated GND density as the deformation becomes sufficiently large. 


2.3. Statistical model of dynamic recrystallization nucleation 

During dynamic recrystallization (DRX), new grains are initiated according to complicated and still 
unclear atomistic mechanisms and usually correlated with stress/strain fields at continuum level (Galiyev 
et ah, 2001; Ghauve et ah, 2015). Together with considering the experimentally revealed grain boundary 
(GB) bulging mechanism (Ponge and Gottstein, 1998; Wusatowska-Sarnek et ah, 2002; Miura et ah, 
2007) , we make the following core assumptions: (a) DRX grains nucleate when some material at a 
GB undergoes certain atomistic transformations driven by thermal fluctuation and local stress to form a 
stable nucleus; (b) the time scale of this transformation and the formation of a stable nucleus is at atomic 
level and hence instantaneous as compared to the applied deformation, and (c) the length scale of this 
transformation is small relative to the FFT-GP computational grid. These consideration are similar to 
those made by Niezgoda et al. (2014) in modeling of the stochastic twin nucleation, and is illustrated in 
Fig. 2(a) to show schematically the length scale separation. 




(a) (b) (c) 

Figure 2: (a) Schematics of the incorporation of GB bulging mechanism in the current DRX nucleation model. 
Heterogeneous microstructure is sampled using a regular FFT-CP grid, which can be divided into many sub-cells 
that can undergo discrete atomistic transformations. Based on our statistical modeling, the local dislocation 
density p(x), defined on FFT-CP grid, is linked with our predicted (b) probability of DRX nucleation and (c) 
the corresponding probability density function of the nucleation strength. 


To formulate a statistical model, we assume that the underlying atomistic transformations, occurring on a 
finely divided subgrid at each FFT-GP gridpoint (Simmons et ah, 2000), are independent and identically 
distributed and that a stable nucleus will form when some minimum number of subgrid points, denoted 
n*, transform. Further we introduce a quantity k{x.) at material point (FFT-GP grid) to measure the 
nucleation propensity or the driving force for transformations, which, as to be discussed below, will be a 
function of the local dislocation content that can be computed from FFT-GP. Given that other factors 
that control nucleation are not explicitly included, we expect that the number of transformation events, 
N, is not strictly deterministic, but rather is a random discrete variable. Our goal is then to link N and 
k{x.) by formulating a continuous distribution, which gives the probability of a stable nucleus forming at 
a specific GB (FFT-GP) gridpoint given /c(x), and then design a stochastic process accordingly for model 
implementation. We leave the detail of derivation in Appendix A and present directly the probability 
of nucleating a DRX grain at a FFT-GP gridpoint (assuming n* = 1) 

P(N(fc) > 1) = (9) 
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where kc 


C~^ exp ^ ^ with C and q being material constants and Qdrx the apparent activa¬ 

tion energy associated with DRX (Ding and Guo, 2001; Rollett et ah, 2004). To further appreciate Eq. 
(9), we introduce the concept of a stochastic nucleation strength k(x), which is a local material property 
and depends on the subgrain/atomic structure. Nucleation will occur when the local propensity exceeds 
the local strength, i.e., k(x) > k(x). The probability density of k can be determined as 


fK(k) 


dP {k < k) _ dP (N(fc) > n*) _ ^ -{^y 

dk dk kr 


( 10 ) 


The inherent stochastic nature we introduce here is the key distinctive feature compared with previous 
works (Ding and Guo, 2001; Takaki et ah, 2008, 2014). 

The construction of the above statistical nucleation model is independent of the exact expression of k{x). 
Here we simply use the total dislocation density 

Af 

fc(x) = /9tot(x) = + /3 sSd(x) +PgNd(x) (H) 

a=l 


as the relevant measure of the DRX nucleation propensity. As the mechanistic details of nucleation are 
clarified, new forms of k{x) can be incorporated without changing the overall framework. Figs. 2(b) and 
2(c) plot, respectively, the probability of DRX nucleation as a function of dislocation density (Eq. (9)) 
and the probability density function (p.d.f.) of the nucleation strength (Eq. (10)), for the simulation to 
be described in Sec. 3. The parameters kc and q were fit to the macroscale stress/strain curve shown 
in Fig. 4(a). Note that as nucleation only occurs at a small fraction of potential sites, kc must be 
significantly larger than the average dislocation density. This relationship is explored in more detail in 
Sec. 4.1. 


To implement the above statistical nucleation model, we initially draw a nucleation strength k(x) for 
each EFT gridpoint x from the Weibull distribution (i.e., Eq. (9)) specified by two parameters q and kc- 
Here n(x) is considered as a static quantity during the deformation until a DRX nucleation event occurs 
at X, which will assign a new value to k(x) due to the change of atomistic environment. The reassignment 


c- 


f Qdrx a 




, where Qdrx 


of K,{x) again requires the specification of q and kc- Recall that kc = 
is the activation energy characterizing the DRX nucleation mechanism and essentially related to the 
nature of the GB. Gonsidering that the initial GB obtained via long term heat treatment is apparently 
quite different from the boundary of recrystallized new grains, we use a different nucleation parameter k^ 
to draw the new k(x). Gurrently it is treated as k^ = Snuci^c with Snuci being simply a fitting parameter. 
In theory, Snuci can be completely determined once the knowledge of Qdrx, <1 and C are available. Since 
experiments have confirmed that in copper the activation energy of migration of low angle GBs is much 
higher than that of high angle ones (Viswanathan and Bauer, 1973), it is expected that s^uci could 
deviate significantly from 1 in reality. Nevertheless a further estimation of s^uci requires knowing the 
absolute value of nucleation rate (related to the coefficient C) and is thus difficult to measure. 


2.4- Phase-field model 

The subsequent evolution of the grain structure containing nucleated DRX grains can be treated by 
extending the PE grain growth model of Ghen and Yang (1994) by including the reduction of the stored 
strain energy, which competes with the reduction of the curvature to determine the migration of the 
recrystallized grain boundaries, as an additional driving force. Moelans et al. (2013) have modified the 
PE grain growth model of Ghen and Yang (1994) by using two order parameters r^rex and r/def to describe 

8 










one grain, with the deformed state corresponding to {ijrex = 0,77def = 1) and the recrystallized state to 
(^rex = l^^def = 0)- Wc further extend the original 2D PF formulation of Moelans et ah (2013) to 3D. 
The total free energy consists of GB energy and the stored strain energy, namely, 

Ftot = -Fgb + i^def = j /gb(x)dx + y/def(x)dx, (12) 

where the local GB energy density /gb (the dependence of x is omitted in what follows) is formulated as 
(Ghen and Yang, 1994; Moelans et ah, 2013) 


Na 


/gb 


= E 


6 wgb 


Vrex.,i I ^def,i ^rex,; 


i=l ^8*’ 


Ng Ng 

4vdei,i I + '^Vdei,j 

U = 1 3>i 


^def,z 3 2 ^^ _ 


2 2 
3 




i=i 


^ g^gb^gb ^(V?7 rex,i f + (Vr?def,i)^) 


(13) 


where Ng is the total number of grain orientations under consideration, cjgb the GB energy that can be 
taken from experiments, and /gb the width of the diffuse interface zone associated with GB in PF method. 
The particular relationships between the coefficients were developed for the quantitative simulation of 
grain growth by Moelans et al. (2008b). The stored energy density /def is formulated as 


/def — -^store (^) 



= /3tot(x)Cli^^ 



(14) 


where the predicted stored energy F/store(x) by FFT-EVP is approximated (Hull and Bacon, 2001) as 
^store(x) ~ ptot(x)CM/>^ where C is a dimensionless constant of 0-1 and the total dislocation density 
is given in Eq. (11). The evolution equations of the order parameters simply follow the Allen-Gahn 
equation (Allen and Gahn, 1979), which in our case takes the following form (Moelans et ah, 2008b): 


^^rex,i = 

Ot 3 y /gb / ^^^rex,i 

^'^def ,i _ f Afgb \ SF tot 

dt 3 V Igh J Sr]def,i ’ 

where Mgb is the GB mobility. To model the property of a polycrystal, the RVE should contain sufficient 
number of grains, which leads to a large number of coupled partial differential equations in Eq. (15). To 
solve this numerical difficulty, we adopt the sparse data structure technique proposed by Gruber et al. 
(2006) and Vedantam and Patnaik (2006) to reduce the effective Ng considered for each PE node. It also 
needs to be pointed out that the PE simulation does not need to reside on the same computational grid as 
in EET-EVP, owing to the different length scales of the corresponding physical processes as indicated in 
Eig. 2(a). A finer PE grid can be used by interpolating the EET-EVP grid to increase spatial resolution 
when simulating the DRX grain growth. To do this, linear interpolation is used to obtain the finely 
sampled smooth fields such as stress, strain, and dislocation densities; the refined grain structure with 
a refinement ratio PE grid/EET-EVP grid = 2:1 is obtained by assigning the orientation value of the 
nearest original grid point to the new interpolating points, which is only involved in grain boundaries 
and junctions. The grain ID after PE simulation is determined by max |{77rex,i, |- 

If a time increment At is used to solve the deformation kinematics (Eq. (1)), the corresponding number of 
PE steps AVpf for microstructure relaxation that can take place before further mechanical loading is sim¬ 
ply given as AVpf = where St is the real time of one PE simulation step, and [xj means the largest 


i = l...N 


9^ 


i = 1. ..Ng, 


(15) 
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integer not greater than x. The determination of St follows Moelans et ah (2008a) where parameters were 
taken so that 6[m^/J]co’gb [J/m^] = /gb [m], /gb = V^Sx, and Mgb [m^/(J • s)] St [s] =0.0375[m^/J]/gb [m] 
where Sx is the discrete grid spacing in PF simulation. Such numerical specification has been shown to 
yield a relative (numerical) error smaller than 5% in terms of simulating GB motion (Moelans et ah, 
2008a). It is then given that St [s] = 0-225|m^/^Jjjyjj/m ] ^ expected is controlled by GB energy 

and mobility. 


In principle, GB energy and mobility should depend on the misorientation and inclination, which leads 
to grain growth in anisotropic systems and has been considered by previous PF models (Kazaryan 
et ah, 2000; Ma et ah, 2004; Suwa et ah, 2007; Moelans et ah, 2008a). Nevertheless, the experiment 
(Wusatowska-Sarnek et ah, 2002) to be simulated and compared shows mostly equiaxed DRX grains, 
unlike what was predicted by PF grain growth in anisotropic system (Kazaryan et ah, 2002); the es¬ 
sentially random texture found in the experiment (Wusatowska-Sarnek et ah, 2002) also frees us from 
considering the high anisotropy in GB mobility that manifests only in highly textured polycrystals (Ma 
et ah, 2004; Suwa et ah, 2007). As a result, we assume isotropic grain growth in the current model and 
take GB energy cjgb and mobility Mgb as constant. 


2.5. Micro structure update and stress redistribution 


Once new grains are formed, state variables, e.g., dislocation density and grain orientation, should be 
updated based on DRX mechanisms. The general picture that DRX grains are “dislocation-free” or with 
a much lower dislocation content describes the dislocation activities occurring at subgrain-level (Sakai and 
Ohashi, 1990) and consequently cannot be directly applied to the dislocation density in the current model, 
which is defined at a level much larger than the subgrains (Fig. 2(a)). The “dislocation-free” nucleus 
at subgrain-level will co-deform with its surrounding grains and will rapidly accumulate a significant 
dislocation density before it grows to the length scale of a typical continuum plasticity simulation. This 
has been confirmed experimentally by transmission electron microscopy (TEM) observations (Sakai and 
Ohashi, 1990; Sakai, 1995). Based on this length scale separation, we update the the GND density of 
material points in new grains as 


PgND ^ '^softenPcND 


(16) 


where 0 < ^soften < 1 is a fitting parameter, while keeping the SSD density unchanged. Physically 
this means that we ignore the very initiation stage of the growth of subgrain nuclei, which has actually 
been implied in our nucleation model (Fig. 2(a)). Numerically this ensures the numerical stability in 
an otherwise high-contrast “mixture” of soft and hard grains. The change of GND accounts for (i) the 
underlying change of subgrains orientation, i.e., dense dislocation cell walls transforming into high angle 
boundaries, and (ii) the elimination of pre-existing GBs. On the other hand, the SSD is expect to be 
more or less the same at the coarse-grained FFT grid. 


Regarding the crystal orientation, experiments have shown formation of mo derate-angle subboundaries 
or twins associated with DRX grains. However, there is no overall correlation between DRX and the 
resulted grain orientations and the texture remains unchanged (or slightly more randomized) during 
DRX (Wusatowska-Sarnek et ah, 2002). For our first demonstrative case-study, we simply assume that 
the crystal orientation of the new DRX grain points will inherit the current local values, leaving the 
influence of DRX grain orientation to future studies. 


Once the nucleation step is complete and the state variables updated, subsequent PF relaxation will drive 
the new grains to consume the surrounding old grains. The change in local dislocation density and lattice 
orientation (the part due to deformation) are accompanied by changes in the elastic strain fields and 
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the stress needs to be re-calculated before the next deformation increment. If a new equilibrium stress 
state is not calculated, the resolved stresses in the now softer recrystallized grains will be significantly 
higher than the slip resistances, resulting in non-physical high strain rates. In order to maintain the 
compatibility, the previous strainrate field (before DRX) is prescribed and new equilibrium stress field 
due to the change of microstructure is calculated before further mechanical loading is applied in FFT- 
EVP. The deformation continues following the dynamically coupled scheme described above until the 
prescribed target strain increment is achieved, as shown in Fig. 1(b). It needs to be pointed out that once 
DRX has been initiated, the nucleation event occurs so frequently that the PF microstructure relaxation 
occurs virtually every FFT-EVP step. As a result, during the implementation of Eig. 1(b), DRX grains 
that have been grown for a few steps are still allowed to grow as long as there is a driving force due to 
the store energy difference. 


3. Results 

Simulation was carried out after the experimental study of Wusatowska-Sarnek et al. (2002). The material 
was 99.99% pure Cu (4N) with equiaxed grains exhibiting an average size of ~ 230 /im. The texture 
was nearly random with weak (< 2 times) random (111) and (10 0) components. The representative 
volume element (RVE) in our simulation was produced by DREAM3D (Groeber and Jackson, 2014). 
The RVE contained 191 grains and the grid size Iq = 21.4 jam was chosen such that the initial average 



(a) (b) (c) (d) 


Figure 3: The input grain structure is randomly colored in (a) with the corresponding inverse pole figure (IPF) 
shown in (b) (the preferred orientation is parallel to the compression axis). IPFs at 50% strain stages are shown 
for (c) a standalone FFT-EVP and (d) an integrated (FFT-EVP + PF) simulations at 723 K, of which the 
stress-strain curves are shown in Fig. 4(a). 


grain size matches the experimental value. The initial grain structure and inverse pole figure (IPF) are 
shown in Figs. 3(a) and 3(b), respectively. Experimentally DRX was not observed in polycrystalline 
4N Cu at 473K at strains < 50% (Wusatowska-Sarnek et ah, 2002). Therefore a standalone EET-EVP 
without PE coupling was used to calibrate the dislocation model. The RVE was subjected to simple 
compression at a strain rate of 1.6 x 10“^/s to match the experimental conditions. The temperature 
dependence of the model was verified against the early part of the stress strain curve for 4N Cu at 723 
K before the onset of dynamic recrystallization (See Eig. 4(a)). The resulting parameters are listed in 
Table 1 and the corresponding simulated stress-strain curve is shown in Eig. 4(a). It needs to be pointed 
out that the initial SSD density is taken as 2.5 x 10^^ m“^ at 473 K and 2.8 x 10^^ m“^ at 723 K. 
This is found necessary in order to account for the significant difference between the initial yield points 
of the experimental stress-strain curves shown in Eig. 4(a). One possible reason is that effects such 
as grain boundary strengthening and/or temperature-dependent activation volume have been ignored 
in the current study. The evolution of simulated dislocation densities are shown in Eig. 4(c), which as 
expected shows that the mobile dislocation population is only a small fraction of that of the immobile 
dislocations. 
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Table 1: Parameters of dislocation-based constitutive model for copper. (fFrost and Ashby (1982); fDavis et al. 
( 2001 ).) 


Symbol 

Value 

Meaning 

Qslip 

3.3 X 10-19 J t 

Activation energy barrier for slip 

^bulk 

3.51 X 10-19 J t 

Activation energy barrier for climb 

lo 

21.4 /im 

Grid size for a 64^ RVE containing 191 grains 

e 

200 

Gonstant for effective length of GND 

Cl 

0.5 

Gonstant for passing stress 

C2 

3.0 

Gonstant for jump width 

C3 

3.0 

Gonstant for obstacle width 

C4 

8.0 X lOi" m-i 

Gonstant for lock forming rate 

C5 

10.0 

Gonstant for athermal annihilation rate 

C6 

3.0 X 1019 m-i 

Gonstant for thermal annihilation rate 

Cl 

7.0 X 10-28 m^s^s- 

Gonstant for dipole forming rate 

C8 

0.24 

Gonstant for nonlinear climb 

Table 2: Material properties and parameters of nucleation and phase-field models for copper. (fMurr 
fVandermeer et al. (1997).) 

Symbol 

Value 

Meaning 

^gb 

0.625 J/m^ t 

Grain boundary energy 

Mgb 

145 mV(MJ . s) t 

Grain boundary mobility 

C 

0.25 

Constant for stored strain energy 

kc 

300 [10^2 m-2] 

Characteristic dislocation density for nucleation 

Q 

4.4 

Exponent for Weibull distribution of nucleation rate 

<^nucl 

0.05 

Accounting for the change of kc 

<5soften 

0.9 

Eraction of GND inherited in DRX grains 


With the calibrated constitutive law in hand, we now integrate FFT-EVP with our nucleation model and 
PF model to implement the complete DRX simulation as depicted in Fig. 1(b). Experimental values were 
used for GB energy and mobility (Murr, 1975; Vandermeer et ah, 1997) while the other nucleation and PE 
parameters were fit to the experimental stress/strain data (see Table 2). The resulting fit reproduces both 
the experimental stress-strain curve and hardening rate (HR) evolution at 723 K with good agreement, 
as shown, respectively, in Eig. 4(a) and Eig. 4(b). In the standalone crystal plasticity simulation, the 
HR monotonically decreases to zero due to the dynamic recovery (i.e., annihilation terms in Eq. (6)). In 
the integrated simulation, the additional softening due to DRX starts to produce a slight deviation on 
the HR curve around 10% strain and gradually macroscopic softening, i.e., a negative slope on the stress- 
strain curve (or equivalently negative HR). The corresponding evolution of average dislocation densities 
are shown in Eig. 4(d), where the average density of SSD exhibits a significant decrease during DRX and 
is a major contributor to the softening. As will be shown later, the softening is mainly attributed to the 
local stress relaxation and redistribution in the neighborhood of newly nucleated and growing grains. It 
should be noted that all model parameters were exclusively fit to the macroscopic stress-strain curves, 
and the resulting evolution of dislocation densities is a key prediction of the coupled simulations. 
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True strain 


True strain 


(c) 


(d) 


Figure 4: Comparison between experimental (Wusatowska-Sarnek et ak, 2002) and simulated (a) compression 
stress-strain curves at 473 K and 723 K and (b) the hardening rate at 723 K, with the corresponding evolution of 
dislocation densities at (c) 473 K based on FFT-CP simulation and (d) 723 K based on the integrated (FFT-EVP 
+ PF) simulation. 


Another key prediction is the full-held microstructural information during DRX. Figs. 5 plots the 
evolution of total number of grains in the RVE, the volume fraction of recrystallized grains and the 
mean grain size during the deformation. The overall sigmoidal form of the recrystallized volume fraction 
agrees with the experiments (Blaz et ah, 1983; Rollett et ah, 2004), and the continuous increase of the 
number of gains suggests grain rehnement, which agrees with the experiment as well. The underlying 
grain structure evolution is shown in Fig. 6 by plotting snapshots of grain maps during the compression. 
It can be seen that DRX grains have been formed at both grain boundaries and junctions. Note that the 
grain structure has changed so significantly that the detailed DRX nucleation and growth processes can 
only be achieved by more quantitative analysis of the simulation result, which will be presented in the 
following Discussion section. The simulation also shows that new grains are mostly equiaxed, which is 
consistent with the corresponding experiment (Wusatowska-Sarnek et ah, 2002). IPFs at the macroscopic 
strain of 50% for the standalone FFT-EVP and integrated simulations at 723 K (Fig. 4(a)) are plotted 
in Figs. 3(c) and 3(d), respectively. The exhibited formation of (1 0 1) fiber texture is consistent with 
the experiment (Wusatowska-Sarnek et ah, 2002). In addition, DRX essentially places no difference in 
terms of the texture, which is also consistent with the experimental conclusion that DRX only causes 
slightly weakening of the (1 0 1) fiber texture (Wusatowska-Sarnek et ah, 2002). 
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Figure 5: Evolution of (a) the total number of grains, (b) volume fraction of recrystallized grains and (c) mean 
grain size during the integrated modeling at 723 K. 



(a) (b) (c) (d) (e) 


Figure 6: Grain map, colored by grain ID, at the strain of (a) 0%, (b) 14.4%, (c) 26.3%, (d) 38.2%, and (e) 50.0% 
during the integrated modeling at 723 K. The DRX grains at different strains can be identified by comparing the 
corresponding microstructures in (b)-(e) to that in (a). 


4. Discussion 

4 . 1 . Critical strain of DRX 

The initiation of DRX is well documented to occur before the strain Sp that corresponds to the peak 
stress (Sakai and Jonas, 1984). This threshold strain 5c, known as the critical strain of DRX, is a very 
fundamental DRX parameter for practical purpose. However, it can only be determined under isothermal, 
constant strain rate conditions, involving a large number of interrupted tests (Manonukul and Dunne, 
1999). It is expected that experimental identification of 5c based on 2D cross-section should be always 
larger than the actual strain corresponding to the very first DRX event, for sufficient volume fraction of 
DRX grains is required to yield statistically observable new grains in a randomly chosen 2D cross-section. 
Our integrated modeling provides an accurate way of predicting 5c by simply examining the first DRX 
nucleation event. Table 3 lists our simulation result together with experimental results performed on 
polycrystalline 3N Cu at a slightly different temperature (700 K) with various strain rates. As expected 
our predicted 5c is significantly smaller than the experimental measurement, which should only serve 
as the upper bound of the actual 5c. Our integrated modeling suggests that in reality the initiation of 
DRX may occur at a much earlier stage than what is revealed by the existing experimental technique 
and mean-field theory. In Fig. 7(a), the dislocation density distributions at various strain levels 
prior to DRX nucleation are plotted. With increasing applied strain the dislocation density distribution 
shifts towards towards higher dislocation densities and broadens. After an applied strain of 21.5% 
(^ 5p), however, the mean of the distribution shifts back towards lower values, indicating a macroscopic 
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Table 3: The critical strain to initiate DRX in polycrystalline copper. Experimental data are taken from 

(Manonukul and Dunne, 1999). (^The simulation prediction should always be smaller than the experimental mea¬ 
surement.) 


Strain rate (s ^) 

Experiments, 700 K 

Simulation, 723 K 

5 X 10“^ 

5 X 10-3 

5 X 10-2 

1.6 X 10-3 


0.15 

0.19 

0.23 

0.061' 

£p 

0.22 

0.30 

0.36 

0.22 


0.68 

0.63 

0.64 

0.27t 



1.6 



CO 

1=5 

a; 


o 

^ • 

n®^_^^^_I 

^ 0.1 0.2 0.3 0.4 0.5 

True strain 


Figure 7: (a) The dislocation density distributions at various strain levels prior DRX events and (d) the corre¬ 
sponding expected nucleus density (i.e., volume percentage of DRX nucleation at mesoscale voxel level). The 
nucleation strength distribution (Fig. 2(c)) is also plotted in (a). 


softening. In addition, the distribution of dislocation density overlaps with the distribution of nucleation 
strength only at tails (Fig. 7(a)), suggesting that it is the material points with extreme values of 
dislocation density that actually controls the initiation of DRX. Integrating the overlap between the 
two distributions, f f{k{x)) • f^{k)dk^ gives the expected nucleus density (i.e., volume percentage of 
mesoscale voxels undertaking DRX nucleation) at the given deformation state as shown in Fig. 7(b), 
which is a direct measurement of instantaneous nucleation rates. It suggests that the DRX nucleation 
rate increases significantly once initiated and reaches a maximum at ^ followed by a gradual decrease. 
As discussed above, DRX is triggered by extreme sites, which keep increase as the material continue to 
work-harden. Initially, the work-hardening outweighs the softening due to recovery and DRX, and the 
nucleation rate reaches the maximum at ^ at which the work-hardening is balanced by the softening 
and the apparent hardening rate becomes zero as shown in Fig. 4(b). The post-Cp softening, originated 
from the stress redistribution that will be discussed in Sec. 4.4, implies that the population of extreme 
sites will no longer increase (as the dislocation density distribution is shifting towards lower values as 
shown in Fig. 7(a)) and a macroscopic stress-drop is developed as the dynamic system approaches the 
steady-state. 
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4-2. Nucleation and growth of DRX grains 


The growth of the very first DRX grain is shown in Fig. 8. This DRX grain is formed at a normal grain 
boundary. In addition, the new grain grows in a wedge-like fashion to maintain a triple line with the two 
old grains, wherein the equilibrium triple junction configuration (120°) does not need to be reached. For 
one thing, grain growth is now mainly driven by the stored energy difference rather than the curvature; 
for another, the dynamic nature of DRX is very likely to prevent the equilibrium configuration from 
being reached before other nucleation and growth events come into play. 






(a) 


(b) 


(c) 


(d) 


Figure 8: Growth of the first nucleated DRX (red grain) during the PF relaxation of the integrated modeling at 
PF time t = (a) 0 s, (b) 0.016 s, (c) 0.021 s, and (d) 0.030 s. Here the initial stage (a) corresponds to the global 
(physical) time 37 s and the time step in FFT-EVP simulation is At = 0.03 s. 


Apart from nucleation at grain boundaries, experiments have indicated that triple junctions may be the 
preferential nucleation sites for DRX grains (Miura et ah, 2005). However, as has also been pointed 
out by Miura et al. (2005), from 3D point of view, quadruple junctions should be more preferred over 
triple junctions in terms of DRX nucleation due to the fact that stress or strain can concentrate more 
easily at quadruple junctions owing to its non-equilibrium nature. In our simulation, we have actually 
identified the nucleation of DRX at a quadruple junction, as shown in Fig. 9. A close inspection of the 
3D structure as shown in Fig. 10 suggests that as in the case of nucleation and growth at a normal GB, 
the DRX grain formed at a quadruple junction also tends to maintain a triple line with each pair of the 
neighboring grains by growing in a wedge-like fashion. 



(a) (b) (c) (d) 


Figure 9: Growth of the a nucleated DRX (red grain) at a quadruple junction during the PF relaxation of the 
integrated modeling at PF time t = (a) 0 s, (b) 0.016 s, (c) 0.021 s, and (d) 0.030 s. Here the initial stage (a) 
corresponds to the global (physical) time 47 s and the time step in FFT-EVP simulation is At = 0.03 s. 


It would also be of great interest to analyze all nucleated DRX grains in our integrated simulations and 
perform statistics on the nucleation at GB, triple junctions and quadruple junctions. However, this is 
beyond the purpose of the current paper and we leave that as a topic of future work. 
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(a) (b) (c) (d) 


Figure 10: Individual triple junctions formed by every three contacting grains at a DRX grain formed at quadruple 
junction as shown in Fig. 9. The viewpoint has been adjusted accordingly for the best illustration. 


4 . 3 . Kinetics of DRX 

For static recrystallization, the kinetics can be analyzed using the Avrami equation: 

X = l-exp[-Bt^], (17) 

where X is the volume fraction recrystallized, t is time, m is a constant called Avrami exponent, and B 
is also a constant. During DRX if the strain rate 5 is constant, the time in Eq. (17) is usually replaced 
by t = {e — Sc)Ie and the kinetics can again be analyzed using the Avrami equation (Roberts et ah, 
1979), which is shown in Fig. 11. The Avrami exponent m = 1.42 from our simulation (Fig. 11(a)) 




Figure 11: Avrami analysis of simulated DRX kinetics: (a) the determination of Avrami exponent m and (b) the 
comparison of between simulated and htted volume fraction recrystallized. 


falls in the range of the experimental values 1.0 ^ 2.4 revealed by the data on commercial purity copper 
(Garcia et ah, 2000). 

The slope in Fig. 11(a) shows some variation, indicating that the Avrami exponent is decreasing when 
deformation becomes large. The nonuniform Avrami exponent has also been observed in static recrystal¬ 
lization, where the Avrami exponent decreases as the recrystallization proceeds. This is attributed to the 
lack of uniformity of the stored energy (Rollett et ah, 1989), which is also the case in DRX. The Avrami 
analysis (Fig. 11(b)) indicates that the completion of DRX is ^ 80%, at which level the steady-state 
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flow has just been reached according to the experimental stress-strain curve of Wusatowska-Sarnek et al. 

( 2002 ). 

The estimated mean grain size T^rex during the compression (Fig. 5(c)) shows grain refinement of 
Drex/Do = 0.5 {Do is the initial mean grain size) up to 50% strain, which is consistent with the com¬ 
mon empirical criterion (Sakai and Jonas, 1984) that a single stress peak DRX usually corresponds to 
2T)rex < Dq. The corresponding experiment data (Wusatowska-Sarnek et ah, 2002) shows D^q^/Dq = 0.2, 
measured after a strain of ^ 1.3 when the steady-state has been reached. It is difficult to “extrapolate” 
current simulation results (Figs. 5(a) and 5(c)) and compare with this experimental result. Nevertheless, 
the initially rapid grain refinement has shown a decrease in the refinement rate. In our future work, the 
replacement of FFT-EVP with a finite-strain framework will eventually allow us to approach the steady 
grain size revealed by many DRX experiments (Sah et ah, 1974; Blaz et ah, 1983; Sakai and Jonas, 1984). 

4.4‘ Softening during DRX 

While the analogy between the formation of “dislocation-free” grains and phase transformation has 
been adopted to intuitively explain the softening of DRX, it has to be noted that DRX is essentially 
related to the formation and migration of grain boundaries, which are extended defects rather than 
phases, and extended defects are strongly driven by stress rather than thermodynamics (Zhao et ah, 
2013, 2014). At the length scale of discrete extended defects, the RVE-average stress is expected to 
exhibit successive drops due to the repeated relaxation of local stresses, which, however, may easily be 
averaged out from further coarse-graining at a larger length scale. Eor instance, while dynamic recovery 
due to dislocation annihilation can certainly lead to softening, the effect is usually balanced by the 
hardening due to dislocation entanglement occurring at the same scale (as shown in Eq. (6)), leaving 
the stress-strain curve approaching monotonically a plateau (Rollett et ah, 2004) without showing any 
macroscopic softening, as shown in the standalone EET-CP simulation in Eigs. 4(a) and 4(b). As a result, 
softening (negative tangent modulus) during DRX must be attributed to dislocation population evolution 
due to stress relaxation at a much larger scale than that of discrete dislocation motions/interactions; 
which is limited to the mean free path of dislocation slip (^ l/yptot) and/or the interaction length of 
dislocation annihilation (such as climb). To this end, the stress field prior to the nucleation of the very 
first DRX grain (Eig. 8) is plotted in Eig. 12(a). Compared with the corresponding grain structure in 
Eig. 3(a), it is obvious that stress concentrations mainly occur at the grain boundaries and junctions, 
supporting the GB bulging mechanism adopted in our simulation. (The dislocation density field has 
a similar distribution.) The change of von Mises stress after the growth of the new grain, shown in 
Eig. 12(b), indicates that the most significant changes (> 5 MPa) are localized within a small region 
that corresponds exactly to the region of the new grain as shown in Eig. 8. (Note that Eigs. 8 and 
12 have the same viewpoint.) In addition, the maximum decrease of local stress in Eig. 12(b) is ~ 9 
MPa, which is at the level of the overall macroscopic stress drop as shown in Eig. 4(a). Similar results 
are obtained when analyzing the stress relaxation upon DRX event at the quadruple junction shown in 
Eig. 9. This suggests that enforcing strain rate compatibility immediately after DRX nucleation can 
effectively redistribute the equilibrium internal stress field to maintain concordance with the evolved grain 
structure. The redistributed internal stress can then influence the subsequent dislocation density (both 
SSD and GND) evolution through the mesoscopic constitutive laws (Eqs. (6) and (7)) and eventually 
lead to the macroscopic softening on the stress-strain curve. The exhibited macroscopic instability is a 
dynamic net result developed by a series of successive discrete DRX events subject to continuous loading 
conditions, which is by no means trivial to perceive based solely on static argument such as instantaneous 
dislocation density. The current interpretation of DRX emphasizes the importance of the time and length 
scale associated with the “dislocation-free” concept in the previous DRX description and suggests that 
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Figure 12: Spatial distribution of (a) von Mises stress prior the nucleation of the hrst DRX grain and (b) the 
resulted change of von Mises stress upon the corresponding phase-held simulation (colorbar unit: MPa). 


internal stress redistribution is responsible for the macroscopic softening. The redistribution of strain 
and stress held during DRX nucleation has actually been indicated recently by Chauve et al. (2015) 
using Digital Image Correlation. It is thus of great interest to carry out more detailed analysis as in Fig. 
12 to yield a statistically more reliable correlation between stress redistribution and DRX; this will be 
done in our future work. 


5. Conclusions 

In this paper an integrated full-held modeling scheme that couples plastic deformation with microstruc¬ 
ture evolution has been established, which provides a general framework for investigating thermome¬ 
chanical processes. The key feature of the integrated modeling scheme is the selection of microstructure 
descriptors that will be used in both constitutive laws of plasticity and the nucleation and growth of 
new phases/grains. As the hrst demonstration, a FFT-based elasto-viscoplastic model and a phase-held 
grain boundary migration model were integrated through the combination of a dislocation-based crystal 
plasticity model and a statistical DRX nucleation model, with the density of dislocations (including 
SSD, GND, and mobile dislocations) being the selected microstructure descriptors to couple plastic and 
microstructural evolution. After calibration solely based on experimental mechanical response of DRX 
in polycrystalline copper, the model delivered fruitful full-held microstructural information, which were 
used to provide quantitative description and analysis of the DRX process. Quantitative agreement with 
experiments in terms of the kinetics and softening during DRX were achieved. 

This effort represents several hrsts in the simulation of thermo mechanical processes. First and foremost, it 
is the hrst successful demonstration of full coupling between crystal plasticity and phase-held simulations 
using a spectral (FFT) framework. Secondly, this work also represents the hrst implementation of a 
dislocation based constitutive theory within the FFT-CP framework with applications to polycrystals, 
further demonstrating the utility of the FFT method as a general micromechanical modeling platform. 
Finally, to the best of the author’s knowledge this work represents the hrst case where softening has been 
handled within the FFT-CP as well. 

Despite the fact that DRX is a well studied phenomenon, our integrated simulation makes several new 
predictions awaiting experimental verihcation. Namely we predict the evolution of dislocation densities 
through the early and intermediate stages of DRX. In addition, we predict that the dislocation content 
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of newly nucleated grains must increase rapidly due to the co-deformation of the new grain with its 
surrounding neighborhood. This suggests that the observed macroscale softening is largely due to the 
effect of local stress redistribution on the dislocation density, particularly a drop in the statistically 
stored dislocation content. We further predict a much earlier onset of DRX than what have been 
reported or predicted previously due to the heterogeneous deformation and stochastic nature of DRX 
nucleation. Finally we show that DRX grains formed at both grain boundaries and junctions (e.g., 
quadruple junctions) tend to grow in a wedge-like fashion to maintain a triple line (not necessarily in 
equilibrium) with the old grains. 

While quantitative agreement of our simulations with experiments have been reached, the current model 
has some limitations with respect to the specific DRX simulation explored here. First of all, the small- 
strain framework of FFT-EVP prevents us from simulating the complete DRX process to compare with 
experiments. Secondly, the special orientation of nucleated DRX grains (e.g., forming the X3 twin 
boundary with the parent grain (Wusatowska-Sarnek et ah, 2002)) has not been considered in the current 
model. Finally, the dislocation-based constitutive laws employed and modified in the current model 
can be further improved towards a more physics-based manner to account for features such as GND 
emission/absorption at grain boundaries or interfaces (Sun et ah, 2000) and GND saturation (Kysar 
et ah, 2010; Oztop et ah, 2013). The future adoption of the finite-strain FFT-GP, would allow a complete 
simulation of DRX and more accurate prediction of the steady state mean grain size. The integration of 
finite deformation kinematics with the phase-field is also crucial for the simulation of other technologically 
relevant thermo mechanical processes. 
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Appendix A. Poisson statistical model for DRX nucleation 

Gonsider a unit area of GB with a nucleation propensity /c(x)^ (for clarity, we will omit the dependence of 
x), then due to thermal and atomic structural fluctuations, the number of transformation events N will be 
a random number and expected to increases with increasing k. Given the assumptions stated above, we 
can say that N(/c) obeys an inhomogeneous Poisson process (Niezgoda et ah, 2014). Gorrespondingly, we 
introduce the cumulative hazard function A{k) that encodes the probability of a subcell transformation. 
The probability of observing N(A:) = n sub-cell nucleation events is thus given by 

P{N{k)=n) = ^ exp [—A(/c)] A(/c)’^ (A-1) 

In order to link Eq. A.l to a probability of nucleation we need to determine the minimum number of 
transformations required to form a stable nucleus, n*, and the form of A{k). Recall that we have not 
specified how the division of a EET gridpoint into sub-cells can be done, and due to lack of a specific 
atomistic mechanisms for DRX nucleation it is unclear what physical events A{k) describes. Glearly there 


^The convention followed in the statistics literature of denoting random variables in bold unfortunately conflicts with 
the usual use of bold typeface to denote tensorial or vector quantities. It should be clear from context whether a particular 
variable is a random scalar or a tensor of vector quantity. 
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is a connection between the two. If a finer sub-grid is chosen, then A{k) is expected to have higher absolute 
values compared to the case of a coarser sub-grid, with, however, the actual reaction rate remained the 
same. This non-uniqueness actually allows us to arbitrarily fix the sub-grid for our convenience such 
that n* = 1 and fit A{k) so that the resulting nucleation probability matches experimental data. The 
probability of nucleation can then be given as 

n* — 1 

P (N(fc) > n*) = 1-P (N(fc) < n*) = I - ^ ’ (^'2) 

n=0 


P(N(/c)>l) = 1 — exp (—A(/c)), 


(A.3) 


which will serve as the numerical criterion to check if a FFT gridpoint will undergo DRX nucleation. 
This is the same as the explicit nucleation model of Simmons et al. (2000) used in phase-field simulation 
of phase transformations, as well as the twin nucleation model of Niezgoda et al. (2014). Following the 
previous studies of Ding and Guo (2001) and Rollett et al. (2004), the cumulative hazard function adopts 
an exponential form 


A{k) = Ck^ exjp 


f Qdrx \ 

V ^bT ) 



(A.4) 


where C and q are material constants, kc = 


C-^ exp 

teristic threshold that the nucleation propensity needs to exceed in order to yield appreciable nucleation 
probabilities. Substituting Eq. (A.4) into Eq. (A.3) leads to Eq. (9) in Sec. 2.3. 




I /cbT J 


and can be interpreted as the charac- 
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